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In this paper, we construct high— order hyperbolic residual distribution schemes for 
general advection diffusion problems on arbitrary triangular grids. We demonstrate that 
the second— order accuracy of the hyperbolic schemes can be greatly improved by requiring 
the scheme to preserve exact quadratic solutions. We also show that the improved second- 
order scheme can be easily extended to third— order by further requiring the exactness for 
cubic solutions. We construct these schemes based on the LDA and the SUPG methodology 
formulated in the framework of the residual— distribution method. For both second— and 
third-order— schemes, we construct a fully implicit solver by the exact residual Jacobian of 
the second— order scheme, and demonstrate rapid convergence of 10—15 iterations to reduce 
the residuals by 10 orders of magnitude. We demonstrate also that these schemes can 
be constructed based on a separate treatment of the advective and diffusive terms, which 
paves the way for the construction of hyperbolic residual— distribution schemes for the 
compressible Navier— Stokes equations. Numerical results show that these schemes produce 
exceptionally accurate and smooth solution gradients on highly skewed and anisotropic 
triangular grids, including curved boundary problems, using linear elements. We also 
present Fourier analysis performed on the constructed linear system and show that an 
under— relaxation parameter is needed for stabilization of Gauss-Seidel relaxation. 

I. Introduction 

In many flow simulations, accurate prediction of solution gradients, such as velocity and temperature 
gradients, are essential for design and analysis purposes as they are directly related to the physical quantities 
of interest: the viscous stresses, the vorticity, and the heat fluxes. However, it is widely accepted that accurate 
and smooth solution gradients cannot be achieved with conventional schemes on fully irregular unstructured 
grids. 1,2 In conventional schemes, the gradients are obtained typically with a lower order of accuracy 
(typically through reconstruction) and they are also subject to numerical oscillations on such grids. The 
resolution of this issue is very important for justifying the use of high-fidelity models in engineering design, 
analysis, and optimization, especially for applications involving complex geometries. The ability to predict 
the gradients on irregular grids is even more critical for grid adaptation, a vital technique for efficient CFD 
calculations in high-order methods, 3 because the grid adaptation almost necessarily introduces irregularity 
in the grid. In fact, current practices in grid adaptation often avoid adaptation in certain regions such as 
within boundary layers where grid irregularity has a severe impact on the solution quality. 4 Therefore, 
numerical schemes that can accurately predict solution gradients on irregular grids need to be developed, so 
that the power of grid adaptation can be fully exploited. 

A method that enables the construction of such schemes is the first-order hyperbolic system method, or 
the hyperbolic method for short, which was proposed in 2007 for diffusion. 5 In the hyperbolic method, which 
is fundamentally different from conventional methods, the solution and the solution gradients are simulta- 
neously computed by solving a hyperbolic system for diffusion. The hyperbolic method was then studied 
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for advection-diffusion in Ref. 6 with Residual- Distribution (RD) schemes. Later, the method was demon- 
strated for the compressible Navier-Stokes equations by a second-order Finite-volume (FV) scheme. ' Since 
then, there have been efforts in developing high-order hyperbolic schemes in the finite-volume method, 8-10 
in the active flux method, 11 and in the RD method 12,13 for unsteady computations. 

In this paper, we focus on the development of hyperbolic RD schemes for two-dimensional problems, 
extending the previous work, 5,6 with several important advances. Improved Accuracy: We propose to 
construct a second-order scheme such that it preserves exact quadratic solutions, which can be accomplished 
by the curvature correction technique. 14 The resulting scheme remains compact, and produces significantly 
improved solution gradients over the previous schemes, which do not preserve exact quadratic solutions. 
Third-Order Accuracy: Extending the improved second-order scheme, we construct a third-order scheme 
that preserves exact cubic solutions. The construction requires quadratic least-squares (LSQ) gradients and 
a higli-order source term discretization. Nonlinear Equation: The improved schemes are extended to a 
nonlinear advection-diffusion equation by the preconditioned conservative formulation introduced in Ref. 7. 
High-Order on Linear Elements: We demonstrate that the third-order scheme does not require curved 
elements for curved boundary problems; it gives more accurate solution and gradients than the second-order 
scheme on the same linear grids. This is a significant advantage, because most high-order methods require 
curved geometries to be represented by high-order curved elements; see Ref. 3. Non-Unified Approach: 
Instead of the fully integrated approach of discretizing the hyperbolic advection-diffusion system as in 
Ref. 6, we discretize the advective and diffusive terms separately. This approach will enable the extension 
to the compressible Navier-Stokes equations for which the eigenstructure of the whole system has not been 
discovered yet. Fully Implicit Solver: We construct a fully implicit solver for both second- and third-order 
schemes. For practical applications, explicit iterations considered in Refs. 5, 6 are not efficient enough, and 
a fully implicit solver is needed. The implicit solver is constructed by the exact residual Jacobian of the 
second-order scheme. It converges in at most 10 iterations to reduce the residual by 10 orders of magnitude 
for both second - and third-order schemes. We demonstrate these features for a series of test problems 
involving fully irregular isotropic aud anisotropic triangular grids and curved boundaries. This work serves 
as a basis for the development of high— order multidimensional hyperbolic RD schemes for more complex 
equations such as the Navier-Stokes equations. The extension of the proposed RD schemes to Navier-Stokes 
equations and problems with shocks and discontinuities will be addressed in subsequent papers. 

The paper is organized as follows. In the next section, we describe the basics of the hyperbolic RD 
scheme: residual evaluations, boundary conditions, and an implicit solver. In Section III, we describe the 
construction of the LDA and the SUPG distribution matrices based on a non-unified approach, in which the 
advective and diffusive terms are treated independently. In Section IV, we discuss the accuracy issue of the 
previous schemes, and propose a guiding principle for constructing improved schemes. A new second-order 
scheme is presented in Section V, followed by the extension to third-order in Section VI. The extension to 
nonlinear advection-diffusion equations is explained in Section VII. Numerical results are then presented 
in Section VIII for both linear and nonlinear problems, followed by concluding remarks in Section IX. For 
completeness, Fourier analysis of the constructed linear system is also presented in Appendix A. 

II. Baseline RD scheme for advection-diffusion 

This baseline scheme is the basis for the development of improved second-order and third-order schemes 
discussed later. For the purpose of our discussion and simplicity, we first describe the details for the linear 
advection-diffusion equation. Extension to the nonlinear equation is discussed in Sec. VII. 

A. Hyperbolic advection-diffusion system 

Consider a two-dimensional advection-diffusion equation: 

d t u + a d x u + b d y u = v(d xx u + d yy u) + s(x,y,u), (1) 

where a and b are constant advection speeds, respectively, in x and y directions, v (> 0) is the diffusion 
coefficient, and s(x, y , u ) is a source term. Following Refs., 6, 12 we rewrite the above equation as a first-order 
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hyperbolic advection-diffusion system: 


d T u + a d x u + b d y u = 

u(d x p + d y q) + s(x, y, u), 

( 2 ) 

II 

{d x u — p)/T r , 

( 3 ) 

H 

( d y u- q)/T r , 

( 4 ) 


where r is a pseudo time, T r > 0 is the relaxation time, and s(x,y,u) is a sum of s(x,y,u) and a physical 
time-derivative term discretized by an implicit time-stepping scheme (see Refs. 12,13 for more details). In 
the vector form, the system is written as 


<9U a <9U b <9U 

dr dx dy 

where 
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( 5 ) 


( 6 ) 


The system is hyperbolic in the pseudo time as shown in Refs. 6. Towards the pseudo steady state, the 
variables p and q approach the solution gradients u x and u y and hence the above equation becomes identical 
to the original advection-diffusion equation with the physical time-derivative discretized by an implicit 
time-stepping scheme (see Ref. 12 ). This is true for any T r , and thus T r is determined not by physical 
constraints but by optimal steady convergence criteria, 6 leading to a non-stiff hyperbolic formulation for 
diffusion. Note that simply dropping the pseudo time derivative will also recover the original equation. In 
this paper, we focus on steady problems, but the resulting steady schemes can be made time accurate by 
including a physical time derivative term in the source term, s( x, y , u), as described in details in Refs. 12- 13 In 
either case, an efficient steady solver is required, and its development is one of the objectives of the present 
paper. 


B. RD scheme: cell residual, distribution matrix, nodal residual 

We discretize the hyperbolic advection-diffusion system on unstructured triangular grids. The domain is 
divided into a set {E} of arbitrary triangular cells (or elements), and an associated set {J} of nodes (or 
vertices). The total number of nodes is denoted by N. We store the solutions (uj. pj, q y ) at each node 
j £ {J}. In the RD method, we first evaluate the residuals over the elements and then distribute the 
residuals to the nodes with a distribution matrix. The cell-residual over an element E £ {E} is defined as 


3 



Figure 1. Schematic of the residual distribution to local nodes and definition of the inward unit normals 
(not-to-scale) . 
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We assume a piecewise linear variation of U (i.e. Co + c x (x — x c ) + c y (y — y c )) over the element, which 
interpolates the three nodal solutions, and perform the integration to obtain 

3 

= - J2 K ‘ u * + Q E dn E , (8) 

2=1 

where i is the local vertex counter (i.e., i = 1,2, 3) for the element E , d£l E is the element area, and 

1 1 3 
Ki = -{Ah^+Bh iy )\ni\ = -A ni \ni\, Q S = |^Q i; . (9) 

i— 1 


Note that n 


i = is the scaled inward normal to the face (edge) opposed to the vertex i, |rq| = 

(see Fig. 1), and n, = (nj x ,rh y ) is the unit inward normal. We note that for every element 


E(=i n < = o and therefore we have E?=i K, = 0. 

The next step is to distribute the residuals to the three vertices by the distribution matrices, Bf , Bf", 
Bf, as illustrated in Fig. 1. The discussion of the distribution matrix is one of the key contributions of the 
present paper, but we leave the choice open at this point. We only mention here that B f is a 3 x 3 matrix, 
which sums up to the identity matrix over the element for conservation. The distribution process results in 


the nodal residual at node j: 


ReSj = 


— y 

3 ££{£,} 


E 




(10) 


where dflj is the median dual volume around the node j (see Fig. 2), and {E :j } is a set of triangular elements 
sharing the node j. The nodal residual is an approximation of the spatial part of the target equations, and 
thus it leads to a semi-discrete form: 

^ = Res.. (11) 

dr 

It can be integrated to the steady state by explicit pseudo-time stepping schemes as in Refs. 5,6, which is 
significantly faster than conventional explicit schemes because the hyperbolic formulation eliminates a typical 
diffusion constraint, 0(h 2 ), where h is a representative mesh spacing, on the explicit time step. 5,6 However, 
it still requires a large number of iterations, especially for anisotropic grids. To improve the convergence, we 
develop an implicit solver in this work, which is explained in the next section. 



Figure 2. Schematic of a median dual volume around the node j. 


C. Implicit solver 

To solve Eq. (11) for the pseudo steady state, we drop the pseudo-time derivative, and define the global 
system of steady residual equations, 


0 = Res, (12) 

which consists of the right hand side of Eq. (11) for all nodes. Note that the resulting residual equation has 
become consistent with the advection-diffusion equation (1) because the pseudo-time derivative has been 
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dropped. We solve the system by Newton’s method: 


U*+! = u‘+AU i , 


(13) 


where U = (u\,p\,q\, U2,P2,Q2, ■ un,Pn,Qn) and l is the iteration counter. The correction AU ! = 
U /+1 - U' is determined as the solution to the linear system: 


<9Res 

dU 


AU* = Res z (U'), 


(14) 


where Res ( is the steady residual vector evaluated by U z . The Jacobian matrix <9Res/$U is exact and 
sparse because the spatial discretization is compact. For each node j e { J}, it involves ( k + 1) 3x3 blocks, 
where k is the number of immediate neighboring nodes to the node j . For example, k = 6 for the node j 
shown in Fig. 2 and therefore, for this particular node, there are seven 3x3 blocks. For a node j, we have 

k 

— JjAUj- - Y, = R^(U'), (15) 

m = 1 


where 


■h = 


dUeSj 




9Resj 

dUm 


(16) 


We may analytically evaluate the diagonal and off-diagonal entries of the Jacobian matrix, i.e., J ; and J j m , 
but this process is rather tedious for general time-dependent advection-diffusion problems, and even more 
difficult for complex systems such as the Navier-Stokes equations. To overcome the difficulty, we implemented 
an Automatic Differentiation (AD) tool based on an operator-overloading technique to evaluate the exact 
Jacobians numerically through chain rule. Thus, the Jacobian matrix is exact up to the round-off error 
for the baseline and second-order schemes. The linear system is relaxed by the sequential Gauss-Seidel 
relaxation with an under-relaxation parameter introduced for stabilization based on Fourier analysis of the 
linear system (see Appendix A for more details on stability analysis) . Typically, the relaxation is performed 
until the linear residual is reduced by five orders of magnitude with a maximum relaxation steps of 1000. 


D. Implicit boundary condition 

Here, we discuss two Boundary Condition (BC) types: Dirichlet (i.e., u is known), and Neumann (i.e., d n u 
is known). We note that because gradients are also the primary variables in the hyperbolic method (i.e., not 
reconstructed from the solution variable), the Neumann type BC is treated similarly as the Dirichlet type 
BC. 



Figure 3. Schematic of boundary nodes for formulation of implicit boundary condition. 


For a Dirichlet type BC, the following set of equations is imposed at the boundary nodes: 


u — Ub 

(p, q) ■ h 

( d x u - p, dyU - q) - fib 


0, 

d s Ub , 
0, 


(17) 

(18) 
(19) 


where Ub is the given boundary value and d s Ub is the derivative of the Ub along the boundary, which can 
be computed with the given Ub . The fib and tb are, respectively, the unit normal and tangent vectors at 
the boundary nodes (see Fig. 3). For a strong Dirichlet type BC, we specify the BC value at the boundary 
nodes and solve the hyperbolic system accordingly. Here, we discretize and solve Eq. (19) according to 
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the hyperbolic RD scheme described in the previous section. That is, the residual vector at the node jb is 
redefined as 


Uj b Ub 


(Ph i Qjb ) ' tb ~ d s u b 


(20) 


Res Jb (2)n fcx +Res j6 (3 )h by 


where Res J6 (2) and Res Jt (3) are the nodal residuals of the second and third equations of the first-order 
hyperbolic system (i.e. , Eq. 3 and Eq. 4, respectively), computed as described in the previous sections. 

For a Neumann type BC, we impose the following equations to the boundary nodes: 


d T u + ad x u + bdyU 

= v(d x p + d y q) + s, 

(21) 

(p, q) ■ n b 

— 

(22) 

) x u - p, d y u - q) -t b 

= o, 

(23) 


where d n Ub is the given gradient of the u on the boundary in the normal direction, and therefore is known. 
For the strong formulation of a Neumann type BC, we discretize and solve Eqs. (21) and (23) according to 
the RD scheme; that is 


Res, i,(l) 


( Pj b ,Qj b ) ’ n b - d n u b 


(24) 


[ Res ib (2)4 s + Res jb (3)t by \ 

where ReSj 6 (l) is the nodal residual of the first equation of the first-order hyperbolic system computed as 
described in the previous sections. 

In either of the Dirichlet or the Neumann type BC, the modified residual is incorporated into the implicit 
solver with the exact Jacobians. In doing so, a care must be taken to avoid zero diagonal entries in the 
Jacobian; the second and third components in the boundary residual need to be exchanged, depending on 
the magnitude of hb x and hb y , to avoid zero diagonals. 


E. Remark on variable advection coefficients 


If the advection coefficients a and b are functions of (x,y), advective fluxes are not necessarily conservative. 
In this case, we estimate the advection vector by the arithmetic average over an element to write 


<9U _ 3U —dU n 

h A b B = Q. 

or ox oy 

where A and B are defined with the averaged advection vector (a, b ) : 


(25) 



a 

v 0 " 


b 

0 -v 1 

A = 

-1/T r 

0 

0 

, B = 

0 

0 

0 


0 

0 

0 


— 1/T r 

0 

0 


The baseline RD scheme is then applicable. If conservative fluxes exist, then the nonlinear formulation 
described in Sec. VII can be employed to construct an RD scheme. 


III. Distribution matrices: non— unified approach 

The RD scheme is defined by the combination of the cell-residuals and the distribution matrix. The 
distribution matrix is mainly responsible for the dissipative behavior and the stability while the cell -residuals 
determine the order of accuracy. Unbounded distribution matrices can also affect the order of accuracy, but 
we consider only bounded ones in this paper. The construction of the distribution matrix often requires 
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the knowledge of the eigenstructure of the target system. The eigenstructure of the hyperbolic advection- 
diffusion system is known 6 and various distribution matrices can be constructed. However, the extension 
to the compressible Navier-Stokes equations is not possible at the moment because the eigenstructure has 
not been successfully worked out yet. In order to develop hyperbolic RD schemes that can be extended to 
the compressible Navier-Stokes equations, we consider a non-unified approach introduced for FV schemes 
in Refs. 7 and 9, respectively, for the compressible Navier-Stokes and advection-diffusion equations. In this 
approach, we construct distribution matrices for advection and diffusion separately and then combine them 
to derive a distribution scheme for the advection-diffusion equation. This approach easily extends to the 
compressible Navier Stokes equations because the eigenstructures of the inviscid terms and the hyperbolized 
viscous terms can be fully analyzed independently. 7 


A. Unified approach 

To illustrate the point of the non-unified approach, we first recall the unified approach of Ref. 6. Consider 
an arbitrary unit vector n = ( n x ,n y ). The unified advection-diffusion flux Jacobian is defined as: 


A n = A n x + Bh y = 


a n ~vn x -vn y 
—ri x /T r 0 0 

— ri y /T r 0 0 


(27) 


where a n = ah x + bn y is the advection velocity projected onto the unit vector n. The Jacobian A n has the 
following real eigenvalues: 


1 


/ „ Av 

, 1 


I „ Av 

2 

a n - < 

V n + Y r _ 

’ Aa "2 

dn + y 

h + Y r \ 


where Ai < 0 and A 2 > 0. The arbitrary but positive relaxation time, T r , is defined as 


T r = 


Ij x 

\J a 2 + b 2 + v/L r ' 


(28) 


(29) 


where the length scale L r is a quantity of 0(1), which may be taken as L r = 1/27T, as used here, or from an 
optimal formula as derived in Ref. 6. The right and left eigenvectors for the Jacobian A n are given by 



— AiT r 

—^2 T r 

0 


-1 /T r 

— ^ 2 ^x 

— A2 fly 

R„ = 

n x 

n x 

-fly 

T _ 1 

’ L ‘ n ~ Ai— A 2 

1 /T r 

Ain x 

Al fly 


fl.y 

fly 

n x 


0 

— (Ai — A2 )n y (Ai — \2)n x \ 


The unified Jacobian A n can be rewritten based on its right and left eigenvectors as 


A n — R-n A n L n , (31) 

where A is the diagonal matrix with the entries of the eigenvalues defined in Eq. (28). The Jacobian matrix 
can be decomposed as (see Ref. 6) 


where 


Ailli + A^ 

,n 2 - A n - 

f A+ 

(32) 

Ai 

Ai A 2 T r 7i x 

AiA 2 T r n y 


fix /Tr 

-A 2 nl 

A2 flxTly 

, (33) 

^"y /Tr 

A2 flx^y 

-A 2 fi 2 y 


~ A2 

- Ai A2 T r fi x 

— \l\2Trfly 


fi x /T r 

Ai n\ 

Al fl x Tly 

(34) 

fly /T r 

\\fl x Tly 

Aih 2 _ 
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Thus, the K; defined in Eq. (9) can be written as 


K t = K- + K+ (35) 

where 

K- = ^Ai ;i n 1;i In^l = ^A- |n;|, (36) 

K+ = l -M, i U 2ti \u i \= l -A+ i \rv i \. (37) 

These matrices are often required for the construction of the distribution matrix, and as can be seen from 
the above, they require the eigenvalues and the eigenvectors, which are not known at present for a hyperbolic 
formulation of the compressible Navier-Stokes equations. 


B. Non— unified approach 


In the non-unified approach, 7,9 we treat the advective and diffusive terms separately. 
Eq. (5) as 

<9U <9F dG „ 

n. 1 " n. 1 " ~K~ ~ Q: 

or ox ay 


Therefore we rewrite 
(38) 


where 



au 


—up 

F = f“ + F d = 

0 

+ 

—u/T r 



0 J 


0 


bu 


-uq 

G = G a + G d = 

0 

+ 

0 


0 


s 

1 

i 


(39) 


(40) 


Similarly, the flux Jacobian A„ is decomposed into two advective and diffusive fluxes, A“ and respec- 
tively: 

A„ = A“+At (41) 

where 



Q"n 

0 

0 ’ 


0 

—ufi x 

— Ufly 

A a — 

0 

0 

0 

A d — 
5 

fox /Tr 

0 

0 


0 

0 

0 


i 

3 > 

0 

0 


(42) 


The distribution matrix is constructed separately for the advective and diffusive terms based on the corre- 
sponding flux Jacobian, and then the two matrices are combined to form a matrix for the advection-diffusion 
equation. This approach can be extended to the compressible Navier-Stokes equations because the eigen- 
structure of the hyperbolized viscous terms is available. 

The advective Jacobian has only one non-zero eigenvalue (A a = a n ) and thus we can simply decompose 
it according to the sign of A a : 


max(0, a n ) 

0 

0 " 


min(a„,0) 

0 

0 ’ 

0 

0 

0 

+ 

0 

0 

0 

0 

0 

0 


0 

0 

0 


(43) 


The diffusive Jacobian has two non-zero eigenvalues, Af = — \/u/T r , \ d = +\fis/T r , with the right and left 
eigenvectors formally in the form of Eq. (30). Note, however, that they are numerically different because the 
eigenvalues are different. The diffusive Jacobian can be split as 


A d T)d \ dj d \ cl TT 

n ~ ± 1 n A n L 'n A 1 A1 1 


A^n 2 = A: 


d + 


(44) 
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( 45 ) 


where III and n 2 are given by Eq. (33) and Eq. (34), respectively, with Ai = Af and A 2 = \ d . 
We define Kf by 


Kf = -Allnd =Kf +K 


where the negative and positive matrices, K) 4 and Kf + are given by 


d+ 


K d 

= 

= Iwl, 

(46) 

K d+ 

2 

= ^A^n 2>i | ni | 

ii 

[O 1 H 

> 

5 

(47) 


Various distribution matrices can be employed for the hyperbolic RD schemes: Lax-Wendroff (LW), 5 
Low- Diffusion-A (LDA), 6 and Streamline-Upwind Petrov-Galerkin (SUPG) distribution matrices. In the 
following sections, we describe the LDA and the SUPG distribution schemes and their formulations in the 
context of hyperbolic RD scheme (the LW formulation is similar to the SUPG). For completeness, both 
unified and non-unified approaches are presented. 

C. Distribution matrix: Low Diffusion— A (LDA) 

The LDA scheme is a multidimensional upwind scheme, meaning that B)' 0 * becomes zero when no signal 
is being sent to the the upwind nodes. In the unified approach, the distribution matrix B f for the LDA 
scheme 15, 16 is defined as 

b; da - k; • (48) 

For the non-unified approach, the distribution matrix should separately be constructed with the K^ a 
and the K + d terms, corresponding to the advection and the diffusion components, respectively. The diffusion 
component of the distribution matrix, Bf, is similar to the Eq. (51), except that the K.f is replaced with 
the K+ d : 

-l 


B“ = K 


+ d 


5>. 

G=i 


+ d 
3 


(49) 


Note that the inverse matrix is usually evaluated numerically, as done here, but the analytical expression for 
B f is also available in Ref. 17, which we used it only for verification purposes. The distribution matrix for the 
advection term, B“, reduces, however, to a scalar equation as given below because there is no contribution 
from the gradient equations: 


B“ = 


P\ 


LDA 0 0 
0 0 


0 0 


A LDA = 


max(0, a Wi )|rtj| 
max(0, a„ .)|? 


(50) 


These two matrices are combined in the following form to define the LDA scheme for the advection-diffusion 
equation: 

B^ DA = (I - I W )B“ + I w Bf , (51) 

where 

uj 0 0 

I,.. = I 0 1 0 , (52) 

0 0 1 

and uj is a weighting function. From a detailed analysis of a one-dimensional hyperbolic advection-diffusion 
system, the weighting function can be evaluated as 6 


uj = 


Re + 2’ 


( 53 ) 
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where we define the Re for the two-dimensional system as 


Re = \/ a 2 + b 2 lv. (54) 

The weights have been introduced for two purposes. One is to guarantee the conservation: 

E B “ A = I > ( 55 ) 

and the other is to ensure that the scheme approaches a pure advection scheme in the limit of Re — » oo, 
and a pure diffusion scheme as Re 0. Note, for example, that /f DA is 0(1) even for a vanishingly small 
advective vector, and thus it will affect the distribution even in the diffusion limit if the weighing function 
is independent of Re. 

In our experience, hyperbolic RD schemes with the LDA distribution matrix, constructed by the non- 
unified approach, have an instability problem in the linear relaxation on some anisotropic triangular grids 
even with the employment of uder-relaxation parameter obtained from Fourier analysis (see Appendix A). 
On the other hand, hyperbolic RD schemes can be successfully stabilized with a suitable under-relaxation 
parameter in the case of the SUPG distribution matrix, which is discussed next. 


D. Distribution matrix: Streamline Upwind Petrov Galerkin (SUPG) 

The finite-element SUPG scheme 18 is known to be applicable to the RD framework. 19 In the unified 
approach, the SUPG distribution matrix for triangular grids is given by 


B? 



(56) 


which consists of the Galerkin part (the first term) and the stabilization part (the second term). In the 
non-unified approach, we construct the SUPG distribution matrix as follows: 


Bf 


= I + D? 


D 


(57) 


where D“ and T)f are the stabilization terms defined independently for the advective and diffusive terms, 


D 


df PG 0 0 
0 0 0 
0 0 0 


df PG 


1 

2 Ef=i m ax(0,a ni )M’ 


D? = 



(58) 


(59) 


Note that the denominator of Eq. (58) cannot vanish unless the advection vector is exactly zero. A small 
numerical value of the order of machine zero may be added to the denominator in order to avoid zero division 
in the case of zero advection vector. Another possibility is to completely remove Df term for zero advection 
vector; both approaches yield identical result. 

It is also important to note that we have Ei=i D“ = Et=i Df = 0, and therefore the conservation 
property is guaranteed: 

E B ? UPG = L ( 6 °) 

i = 1 

A weighting function can be introduced in the stabilization terms of Eq. (57), 

Bf PG = b + (1 _ w)D? + wDf, (61) 
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where co is defined such that w — > 0 as Re — > oo , and w — >■ 1 as Re — > 0, so that the scheme will properly reduce 
to pure advection and diffusion schemes as Re — » oo and Re — )• 0, respectively. The weighting function u> can 
be evaluated, based on a detailed analysis for a one-dimensional hyperbolic advection-diffusion system, 6 as 


oj = 


2 

Re + 2’ 


where we define the Re for the two-dimensional system as 


Re = a 2 + b 2 / v. 


(62) 


(63) 


In our study, however, we observed no noticeable change in the results between Eq. (57) and Eq. (61) with the 
above weighting function. Therefore, the results shown in this paper are all based on the SUPG distribution 
matrix without the weighting function. Details of the effects of the weighting function is a subject for 
future study. We remark that unlike in the LDA formulation, where the weighting function is required to 
preserve conservation, the SUPG formulation does not require any weighing function because conservation 
is guaranteed regardless. 

As can be seen from above, the SUPG distribution matrix is constructed based on a decomposition of 
the stabilization term into advective and diffusive stabilization terms. Therefore, the SUPG distribution 
matrix can be extended to a hyperbolic formulation of the compressible Navier-Stokes equations, for which 
the inviscid and viscous terms can be analyzed independently. 7 


IV. New design principle for advection-diffusion equation 

An important feature of the RD schemes critical to the accuracy condition is the ability to preserve exact 
polynomial solutions on arbitrary grids. Historically, RD schemes have been designed to preserve linear exact 
solutions for hyperbolic systems of conservation laws. This is called the linearity preservation. Such schemes 
are known to yield second-order discretization errors. The baseline hyperbolic RD scheme 5 is constructed 
based on the linearity preserving property. Although it gives second-order accuracy for all variables on 
regular triangular grids, 5 it was found later that the order of accuracy of the gradient variables (i.e. , p and 
q) deteriorates to first-order on irregular grids. In a subsequent work presented in Ref. 6, a better scheme 
was proposed, which attempted to improve the accuracy of the gradients by upgrading the evaluation of 
the second and third components of the residual. In this paper, the scheme of Ref. 6 is referred to as the 
RD-JCP2010 scheme. Results given in Ref. 6 showed that the RD- JCP2010 scheme, relative to the baseline 
scheme, improves the order of accuracy of the gradient variables, but no results were presented for the quality 
of the predicted gradients. As will be shown later in this paper, we found that for a different test problem 
and a more general exact solution, the same RD-JCP2010 scheme does not yield second-order accuracy 
for the gradient variables (i.e., p and q ), and also generates severe oscillations in these variables within the 
domain. It appears that the discretization of the source terms (including those arising from the hyperbolic 
formulation) was not fully compatible with the flux terms, and it is not clear to us, at the moment, how to 
improve the source term discretization to develop a fully second-order scheme. To overcome this difficult 
problem, we instead propose to improve the flux terms such that the residual vanishes for quadratic solutions 
for the advection-diffusion equation. The exactness, as will be demonstrated later, has a significant impact 
on the accuracy and the quality of the gradient variables on irregular grids. The construction of such schemes 
is very difficult in general because second-order gradients need to be recovered from second-order accurate 
solutions on irregular grids for the diffusive term. However, the construction is quite straightforward in the 
hyperbolic method since the gradients are computed simultaneously with the solution variable u. The scheme 
also extends easily to third-order by imposing the exactness for cubic solutions. Note that the exactness is 
not a necessary condition for accuracy. For example, the third-order RD scheme proposed in Ref. 20 does not 
preserve cubic solutions, but was shown to produce third-order accuracy on smooth unstructured triangular 
grids. 

In the next two sections, we describe how to construct second - and third-order accurate schemes that, 
respectively, preserve quadratic and cubic solutions on arbitrary triangular grids. 
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V. Improved second— order scheme: RD— CC2 


Consider the cell -residuals for the baseline hyperbolic RD scheme: 


®u=J {~au x - bu v + up x + vq y + s)dxdy = \ ^ (-a ni Ui\n,\ + upifi^ + vq^n^ ) + s E dO. E 


= 7fT I {Ux ~ p)dxdy = 7fT f \ Uin i* - P E<mE ) , 


T, 


r J E 


T r \ 2 


»= l 


= - 

9 T r 


[ (u y - q)dxdy = E ( * V 

J E - =1 


(64) 

(65) 

( 66 ) 


These residuals are constructed based on the assumption that all variables (i.e., u, p, q) are linear over the 
element, and therefore they do not preserve quadratic solutions. Note that quadratic solutions imply linear 
gradients (i.e., u x , u y , p, q), and the above cell residuals are already exact for linear gradients. Therefore, 
to preserve quadratic solutions, the integrals of u x and u y terms need to be exact for quadratic solutions. 
A simple way to accomplish this is to reconstruct a quadratic element by interpolating the solution at the 
midpoint of each side in the element (see Fig. 4). The midpoint value u rn . is estimated by the Hermite 


3 



Figure 4. Reconstructed P 2 element with virtual midpoints. 


interpolation along the side: 


u mi = Ui - - (ApiAxi + AqtAyi) , 

O 


(67) 


where Tii is the arithmetic average of the solution values at the two end nodes, and A( ), ( denotes the difference 
of the nodal values taken counterclockwise along the edge opposite to the node i, e.g., Ax 3 = x\ — 22 - Note 
that the midpoint value u mi is exact for quadratic solutions (and linear gradients). Once the quadratic 
element is reconstructed, the cell-residuals are evaluated as line integrals with Simpson’s rule applied along 
each side. This reconstruction approach was first introduced in Refs. 21,22, and later employed in the form 
of the Green-Gauss gradient with a high-order curvature correction term . 6, 14,23 In this work, we employ 
the latter implementation to upgrade the baseline residuals. Note that the formula (67) usually requires 
LSQ gradient reconstruction to obtain the nodal gradients , 14, 21-23 but it is not necessary in the hyperbolic 
method because p(= u x ) and q(= u y ) are already available at nodes as a part of the primary variables. 

Applying the curvature correction, we modify the cell-residuals of the baseline RD scheme as 


1 3 v 3 

= ~7; 1 ^2a ni (ui + S^)\n i \ + ^^2(p i n i!Ii +qi n i y ) +s E dfl E , 


i = 1 


i—1 


= 


Q 


7 fr ( \ + ^)n ix -p E dn E ) , 


T r \ 2 

-f 1 

Tr 2 


i = 1 
3 


^2 (Ui + Sf)n iy - q E dn i 


i= 1 


( 68 ) 

(69) 

(70) 
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(71) 


where the curvature correction term, Sf, is defined as 

Si = \ ( A pi Axi + AqiAy z ). 

o 

Note that in Ref. 6, 5“ was mistakenly shown with a negative sign. One may find it convenient that the 
compact RD-CC2 scheme can be implemented as an add-on to the baseline hyperbolic RD scheme: 


<J> 


E 






— ~ 
T r i 2 


•I s = 


1/1 

H 

9 Tr 2 



(72) 

3 \ 


> 

(73) 

<= 1 / 


3 \ 



(74) 


The modified cell-residuals are distributed by the SUPG distribution matrix. The RD scheme based on 
these cell residuals is referred to as the RD-CC2 scheme. This scheme is compact because no LSQ gradient 
reconstruction is necessary. The exactness for quadratic solutions has been confirmed numerically, and are 
compared with the baseline RD and RD-JCP2010 schemes. These results are summarized in Table 1, 


Table 1. Exactness of the hyperbolic RD equations for exact quadratic and cubic polynomial solutions on 
irregular grids. The last column indicates whether the exactness (/) is based on cell and/or nodal residuals. 


Hyperbolic RD Scheme 

Quadratic solution 
u eq. p eq. q eq. 

Cubic solution 
u eq. p eq. q eq. 

Cell/Nodal 

Baseline 

X 

X 

X 

X 

X 

X 

none 

RD-JCP2010 

X 

/ 

/ 

X 

X 

X 

cell 

RD-CC2 

/ 

/ 

/ 

X 

X 

X 

both 


The truncation error orders for the RD-CC2 scheme, the baseline RD scheme, and the RD-JCP2010 
scheme are provided in Table 2. The truncation error orders have been determined by substituting a smooth 
exact solution (taken from Ref. 24) into the residuals for a series of irregular triangular grids. Clearly, 
all schemes have the same truncation error order for the nodal residual. A known theory for hyperbolic 
conservation laws 25,26 states that the discretization error of 0(h 2 ) leads to the truncation error of 0(h). 
However, our numerical results indicate that the converse is not necessarily true. As we will show later, only 
the RD-CC2 scheme, which has 0(h ) truncation order, gives discretization error of 0(h 2 ). 

Table 2. Comparison between numerical truncation error of the nodal residuals (i.e., T.E. = JT B b # b ) 
obtained with RD— CC2 and non— genuine second— order RD schemes on irregular grids. 


Hyperbolic RD Scheme 

Nodal Residuals 

u equation p equation q equation 

Baseline 

0(h) 

0(h) 

0(h) 

RD-JCP2010 

0(h) 

0(h) 

0(h) 

RD-CC2 

0(h) 

0(h) 

0(h) 


We remark that the RD-CC2 scheme is different from the RD-JCP2010 scheme, 6 which has the curvature 
correction only in and residuals. Thus, the whole residual vector of RD-JCP2010 is not exact for 
quadratic solutions as confirmed in Table 1. The RD-CC2 scheme, on the other hand, has the curvature 
correction applied to all the residuals, and therefore is exact for quadratic solutions (and linear gradients). 

To solve the discrete equations, we employ the implicit solver described in Sec. C. The RD-CC2 scheme 
is compact and therefore we linearize the residual vector exactly by using the AD tool. The implicit solver 
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is thus Newton’s method and converges to machine zero (i.e. , ten orders of magnitude reduction in all the 
residuals) typically within 3-5 Newton iterations. 


VI. Third— order scheme: RD— CC3 


To extend the compact RD-CC2 scheme to a third-order scheme, we modify the cell-residuals such that 
they preserve exact cubic solutions and quadratic gradients. Following the procedure explained in Sec. V for 
the RD-CC2, and thanks to the fact that both the Hermite interpolation and Simpson’s rule are exact for 
cubic functions, we employ, again, the curvature correction technique. Extension to third-order, therefore, 
requires two additional steps: 1) high-order source term discretization, and 2) quadratic LSQ gradients for 
p , q , and the source term s. The high-order discretization of the source term is a critical piece for achieving 
third-order accuracy. We derive a formula by estimating the midpoint value by the Hermite interpolation 
and then exactly integrating the source terms assuming a quadratic variation over the element: 


where 


I E 


s dxdy 
p dxdy 
q dxdy 


dnE ’ 

dn E , 

dnE . 


Sf = ^(As Xi Axi + As yi Ay z ), 
Sf = * ( Ap Xi A xi + Ap y . A yi), 

S q = ^ ( Aq x . A Xi + Aq y . Ayi ) . 


(75) 

(76) 

(77) 

(78) 

(79) 

(80) 


Using these high-order source discretizations, we construct the cell -residuals that preserve cubic solutions 
(and quadratic gradients) as follows: 


^ ^ E ( a "A“H - v{sy Xi + s q n Vt ) - -stdn E ) , 


*p = ^ + + ^ dnE 


2 T, 

(j)^ = 4 — 

q q 2 Tr 


1=1 

3 


E U 


i = 1 


’> yi + 


which can also be written in the form of an extension to the compact RD-CC2 scheme: 

3 


E 


<F 




<I> E 

Q 


= \ J2 ( U ( 6 i n Xi + 8 i n Vi) + d SSidflE ) ’ 


1 




4T r 

$ E + — 
q 4 T r 


i = 1 
3 




i = 1 


(81) 

(82) 

(83) 

(84) 

(85) 

(86) 


The curvature correction terms require quadratic LSQ gradients for p, q, and s, but not for u because 
p(= u x ) and q(= u y ) are already available at nodes. Again, the cell-residuals are distributed by the SUPG 


14 of 40 


American Institute of Aeronautics and Astronautics 



distribution matrix. The hyperbolic RD scheme based on these cell residuals is referred to as the RD-CC3 
scheme. The exactness for cubic solutions has been confirmed numerically as shown in Table 3. The results 
of a numerical truncation error analysis for the RD-CC3 are provided in Table 4 and compared with the 
truncation error of RD-CC2. 

Table 3. Exactness of the hyperbolic RD equations for exact quadratic and cubic polynomial solutions on 
irregular grids. The last column indicates whether the exactness (/) is based on cell and/or nodal residuals. 


rr i ,• r< , Quadratic solution Cubic solution „. 

Hyperbolic RD Scheme Uiemental/Nodai 

u eq. p eq. q eq. u eq. p eq. q eq. 

RD-GC2 ///XXX both 

RD-CC3 ////// both 


Table 4. Comparison between numerical truncation error of the nodal residuals (i.e., T.E. = JV B b # b ) 
obtained with the RD—CC2 and RD— CC3 schemes on irregular grids. 


Hyperbolic RD Scheme 

RD-CC2 

RD-CC3 


Nodal Residuals 

u equation p equation q equation 
0(h) 0{h ) 0{h) 

0(h 2 ) 0{h 2 ) 0(h 2 ) 


The discrete equations are solved by the implicit solver described in Sec. C. We do not attempt to 
linearize the RD-CC3 scheme exactly, but simply use the exact Jacobian of the compact RD-CC2 scheme. 
The implicit solver is therefore not precisely Newton’s method, but it converges so rapidly that fully converged 
numerical solutions (i.e., ten orders of magnitude reductions for all residuals) are obtained typically within 
10-15 Newton iterations as will be demonstrated later in Sec. VIII. 

Remark: The RD-CC2 and RD-CC3 schemes with the SUPG distribution matrix can be considered as 
economical and powerful alternatives to higlr-order finite-element methods. 2 ' They are economical mainly 
for three reasons: 1) the number of linear relaxations in the implicit solver increasers linearly, not quadrati- 
cally as typical for diffusion problems, with grid refinement. 2) the second-order advection-diffusion scheme 
is compact, and thus allows an efficient construction of Newton’s method. Furthermore, the same compact 
Jacobian still yields rapid convergence (comparable to Newton’s method) for the third-order scheme. 3) 
second- and third-order accuracy can be achieved (as will be demonstrated later in Sec. VIII) on linear 
triangular elements without introducing curved elements. Thus, high-order curved grids are not needed, 
and these schemes are directly applicable to existing grids composed of linear triangular elements. These 
schemes are also powerful in that they are capable of producing highly accurate and smooth gradients, to 
the same order of accuracy as that of the main solution variable u , on isotropic and anisotropic irregular 
grids. These advantages will be demonstrated in Sec. VIII. 

VII. Nonlinear advection-diffusion equation 

In this section, we describe the discretization of a nonlinear hyperbolic advection-diffusion equation. As 
an example, we discuss the construction of the RD-CC3 scheme with the SUPG distribution matrix. Note 
that the compact RD-CC2 is a subset of RD-CC3 scheme and therefore the discussion here also applies to 
the RD-CC2 scheme. A similar procedure can be applied to other distribution schemes. Throughout the 
discussion, we only consider the non-unified approach, which is applicable to more complex systems such as 
the compressible Navier-Stokes equations. 

A. Nonlinear hyperbolic advection-diffusion equation 

Consider the following two-dimensional nonlinear advection-diffusion equation: 

d t u + d x f + d y g = d x (vd x u) + d y (vd y u) + s(x , y , u), (87) 
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where / and g are nonlinear functions of u, and v = v(u). The advection speeds in x and y directions are 
therefore a(u) = df /du and b(u) = dg/du , respectively. Using the preconditioned formulation proposed for 
nonlinear systems in Ref. 7, we construct a hyperbolic system as 


d T u + d x f + d y g 

= d x p + d y q + s, 

(88) 

Tf r 

OrP 

V 

= d x u - p/v, 

(89) 

T 

± r 0 

— O r q 

V 

= d y u — q/ v, 

(90) 


where the variables p and q are, in the pseudo steady state, equivalent to the diffusive fluxes in x and y 
directions, respectively. As before, the physical time derivative can be incorporated into the source term 
s as done in Ref. 13, but we focus on the steady equation here. Note that the system reduces to the 
target equation, i.e. , Eq. (87), in the pseudo steady state. The system is written in the vector form as a 
preconditioned conservative equation, with the preconditioning matrix P: 


where 


„ i0U <9F <9G 

P 1 1 

dr dx dy 


Q, 


' 1 

0 

0 


S 

0 

T r /v(u ) 

0 

, Q = 

-p/v{u) 

0 

0 

T r /v{u) 


_ -q/v(u) _ 


p _ pa _|_ pd = 

/ 

0 

+ 

-p 

—u 

, G = G a + G d = 

9 

0 

+ 

1 

° 

1 


0 


0 


0 


. - u \ 


(91) 


(92) 


(93) 


The flux Jacobians can be obtained by multiplying both sides of Eq. (91) by the preconditioned matrix 
P, and arrive at 


where 


PF, = P — U* = AU a „ p Gy = P^U S = BU y , 



a(u ) 

-1 0 " 


b(u) 

0 -1] 

A = 

-v(u)/T r 

0 

0 

, B = 

0 

0 

0 


0 

0 

0 _ 


_ —v{u)/T r 

0 

0 


(94) 


(95) 


Note that these preconditioned Jacobian matrices (i.e., A and B) are slightly different from the ones obtained 
for the linear hyperbolic advection-diffusion system. Hence, their eigen-structures are also different from 
those given for the linear system in Sec. II. The differences in the eigen-structures influence the formulation 
of the distribution matrices as described next. Below, we first describe the baseline RD scheme for the 
preconditioned nonlinear system, and then upgrade it to higher order by the correction approach. 


B. Baseline RD scheme 


The cell -residual for the preconditioned system is evaluated in two steps. First, we evaluate the cell -residuals 
of the nonlinear equation by integrating right hand side of Eq. (91) over an element E £ {E} to get the 
unpreconditioned cell-residuals tU: 




E 


E 


■q/ E 

J 



G y + Q)dxdy, 



i—1 


G i n iy )\n i \+Q E dQ E , 


(96) 
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The preconditioned cell-residual, which is distributed to the element vertices, is then defined as 


= 


^ U 

p 


= PT» 


E 


( 97 ) 


where P is evaluated by the arithmetic average of the solution U within the element. The matrix K j 
corresponding to the preconditioned system is defined by 


K, = -P 
2 


1^/<9F 


dU 


dG. 

au ? 


n, ; 


= K“ + K? 


(98) 


where the flux Jacobians are evaluated by the arithmetic average of the solution, and K“ and K f are, 
respectively, the contributions from the advection and diffusion components, defined as 


K a = 

1 

1- , 

f dF a 

dG a 

2 P ( 

v au n ^ 1 

dV n 'y 

Kf = 

1- 1 

( dF d 

dG d 

P 

— — n, + 

— — ru 


2 \ 

1 9U 

QU l y 


= ^ A nJ n i|l 


= 2 A nJ n il> 


where 



a ni 

0 

0 ’ 


0 

-hi. 

-h iy 

A a — 

0 

0 

0 

> 

3 ft. 

II 

-Vh ix /T r 

0 

0 


0 

0 

0 


-Vh iy /T r 

0 

0 


and v is the arithmetic average of the diffusion coefficient over the element, and 

df „ dg 


Tin — Tl x 
OU 


du n ' v ' 


(99) 

(100) 

( 101 ) 

(102) 


It is clear that the eigenvalues of the advective and diffusive components of the nonlinear hyperbolic 
ad vect ion-diffusion system are equivalent to the linear system: 


A a = 


\d _ 

'U 1 


A^ = 


(103) 


Similarly to the linear advection-diffusion system, we construct the SUPG distribution matrix of the 
nonlinear system over the element E as: 


B fUPG = Ii + D ? + Df , 


where 


D“ = 


and 


d? UPG 0 0 

0 0 0 
0 0 0 

D? = - 


K' 


2 E?=iK? + 


2 E?=i max(0, a ni )\ni\ 


, Kf = -A^|n,|, 


(104) 


(105) 


(106) 


where K,f + is constructed based on the projection of the A d onto the A([ running wave. We note here that 
the diffusive flux Jacobian matrix A d has the following eigenvectors: 


( 107 ) 


R-n = 

" — 1/Ai 

— 1/ A 2 

0 

T — 1 

A 1 A 2 
— A 1 A 2 

Al fix 

A 2 fl x 

Al hy 

—\ 2 fly 

Tlx 

Tlx 

— fly 

5 ~ A1-A2 


_ fly 

fly 

h'x 


0 

— (Ai — A2)n y 

1 

(M 

1 
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which is different from the linear case. Thus, the projection matrix II 2 for the nonlinear case becomes 


n 2 


Ai — A 2 


Ai 

-A1A2 fix 
-A1A2 fly 


-A 2 nl -A 2 fi x fiy 

-A 2 n x n y — A2 fi: 


v J 


(108) 


and the K^ + can be expressed as 

K? + = (109) 

The residual at a node j is then defined in the form of Eq. (10). This completes the construction of the 
baseline RD scheme for the nonlinear advection-diffusion equation. 


C. RD CC2 and RD-CC3 schemes 

The high-order schemes, RD-CC2 and RD-CC3, can be constructed by improving the cell -residuals with 
the high-order curvature correction approach discussed in Sec. V and VI. A special attention should be 
made, however, in constructing the high-order hyperbolic RD schemes for nonlinear equations: the A p and 
A q terms in Eq. (71) require the solution gradients, while p and q variables in the nonlinear formulation are 
the diffusive fluxes. With this in mind, the construction is straightforward. 

The cell-residuals of the RD-CC3 scheme for the nonlinear advection-diffusion equation are given by 


= 


- \ E - *0 n *< + w + l s ° dnE 1 


i— 1 


3 

V 


= $ P+2 y T,( 6 >^ + ^ dnE )’ 


i = 1 
3 




%i! v 


where 8f, 5? and <5? are given by Eciuations (78), (79) and (80), respectively, and 


^ (Af Xi Axi + Af Vi Ayi ) , 

Sf = ^ {Ag Xi Axi + Ag Vi Ayi ) , 

< 5 “ = ^ {A(p/u),Axi + A{q/v)iAyi ) , 

= ^{A{p/v) Xi Axi + A(p/v) yi Ayi ) , 
8i /v = ^{A{q/v) Xi Axi + A(q/u) yi Ayi) . 


( 110 ) 

( 111 ) 

( 112 ) 


(113) 

(114) 

(115) 

(116) 
(117) 


For the flux gradients in the curvature correction terms for / and g , we employ the chain rule to avoid the 
LSQ gradient reconstruction: 


f df 3 / 
h = du {rM ’ 


f.f 


f y = q^v = nrM/ v )> 


d A 

du 


(118) 


and similarly for the flux g. The gradients of s, p, q, p/v, q/v, and the source term, s, are computed 
by the quadratic LSQ fit. The RD-CC3 scheme is obtained by distributing the above residuals by the 
SUPG distribution matrix. The RD-CC2 scheme is obtained simply by removing 8f. S’-, Sf, 8 P J V , and 
8 q J v terms from the cell-residuals. Consequently, the RD-CC2 scheme does not require any LSQ gradient 
reconstruction, and therefore it is compact as in the linear case. We confirmed by numerical experiments 
that the RD-CC2 and RD-CC3 schemes preserve exact quadratic and cubic solutions, respectively, by the 
method of manufactured solutions. 
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VIII. Results 


In this section, we verify and examine the accuracy of the hyperbolic RD schemes. In all the test problems, 
we specify the exact solution u on the boundaries, and initialize all the variables in the interior of the domain 
with a zero value, and solve the hyperbolic system for u , p , and q. The linear relaxation is performed with a 
Gauss-Seidel (GS) algorithm to reduce the linear residual by five orders of magnitude with the maximum of 
1000 relaxations. Typically, it takes 300-600 relaxations with the under-relaxation parameter in the range 
of 0.5 to 0.8. The implicit solver is taken to be converged at ten orders of magnitude residual reduction 
for all equations; this is obtained typically with 3-5 Newton iterations for the RD-CC2 scheme, and 10-15 
Newton iterations for the RD-CC3 scheme. 

Simulations were performed using both the unified and the non-unified approaches discussed in Sec. III. 
The convergence speed and the order of accuracy results were found to be nearly identical with either of the 
two approaches. Therefore, we only present the results obtained with the non-unified approach. 

We also note that for some anisotropic, stretched grids, we could not obtain solution with LDA because 
the GS linear solver was unstable. A more efficient and stable linear solver is needed, particularly for the 
LDA scheme, but we did not pursue this. Therefore, we only focus on the hyperbolic RD schemes with 
SUPG distribution method. 

A. Linear advection diffusion problem 

For the linear advection-diffusion equation (1), we verify the order of accuracy and the quality of the 
predicted solution and solution gradients for isotropic and anisotropic irregular triangular grids, and for a 
problem with curved boundaries. 

1. Isotropic Grids 

Consider a steady two-dimensional linear advection-diffusion equation, Eq. (1), with s = 0. This equation 
has an exact exponential solution 28 of the form: 


u(x, y) = C cos(Airri) exp 


/ 1 — Vl + 4A 2 7t 2 jA 
^ 2v 



(119) 


where A and C are arbitrary constants, and £ = ax + by and 77 = bx — ay. This solution is smooth and has 
no singularly on the boundariesand therefore it is appropriate for the order of accuracy verification. 

A series of isotropic irregular grids was generated and the accuracy of the proposed high-order schemes is 
verified against the baseline RD and RD-JCP2010 schemes. The first three coarse grids are shown in Fig. 5. 
Error convergence results are shown in Fig. 6 along with a more detailed information in Tables 5-7. For 
the variable u, all second-order schemes give truly second-order accuracy, and the RD-CC3 scheme gives 
fourth-order accuracy. For the gradient variables, p and q, on the other hand, the RD-CC2 and RD-CC3 
schemes yield, respectively, second- and third-order accuracy, but the baseline RD and RD-JCP2010 do not 
give truly second-order accuracy. 


Table 5. LI error convergence of u for hyperbolic RD schemes with the SUPG distribution on irregular grids 

(A = 2, C = -1, a = 2, b = 1, v = 0.01). 


Grids 

RD 

Order 

RD-JCP2010 

Order 

RD-CC2 

Order 

RD-CC3 

Order 

32 x 32 

5.66E-05 

- 

6.98E-05 

- 

1.91E-05 

- 

2.31E-06 

- 

64 x 64 

1.08E-05 

2.39 

1.31E-05 

2.41 

4.11E-06 

2.22 

1.43E-07 

4.01 

128 x 128 

2.47E-06 

2.13 

2.82E-06 

2.22 

9.87E-07 

2.06 

1.02E-08 

3.81 

256 x 256 

5.73E-07 

2.11 

6.18E-07 

2.19 

2.38E-07 

2.02 

6.94E-10 

3.88 


Iterative convergence results are provided in Table 8 with the total number of Newton iterations and the 
number of GS linear relaxations per Newton iteration. As expected, the implicit solver, which is Newton’s 
method for the baseline, the RD-JCP2010 and the RD-CC2 schemes, converges with 3-5 Newton iterations. 
For the RD-CC3 scheme, the implicit solver, which is not truly Newton’s method, converges also very 
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(a) 8 x 8 


(b) 16 x 16 


(c) 32 x 32 


Figure 5. Samples of the irregular grids used in this study. 



(a) u (b ) p = u x (c) q = u y 


Figure 6. Order of accuracy comparisons of the baseline RD, RD-JCP2010, RD— CC2, and RD— CC3 schemes. 


Table 6. LI error convergence of p (= u x ) for hyperbolic RD schemes with the SUPG distribution on irregular 
grids (A = 2, C — —1, a = 2, b = 1, v = 0.01). 


Grids 

RD 

Order 

RD-JCP2010 

Order 

RD-CC2 

Order 

RD-CC3 

Order 

32 x 32 

5.19E-03 

- 

7.87E-03 

- 

5.25E-04 

- 

7.43E-05 

- 

64 x 64 

1.70E-03 

1.61 

2.04E-03 

1.94 

1.20E-04 

2.13 

5.88E-06 

3.66 

128 x 128 

5.88E-04 

1.53 

5.54E-04 

1.88 

3.18E-05 

1.92 

5.85E-07 

3.33 

256 x 256 

1.99E-04 

1.56 

1.69E-04 

1.71 

9.14E-06 

1.80 

7.06E-08 

3.05 


Table 7. LI error convergence of q (= u y ) for hyperbolic RD schemes with the SUPG distribution on irregular 
grids (A = 2, C = —1, a = 2, b = 1, v = 0.01). 


Grids 

RD 

Order 

RD-JCP2010 

Order 

RD-CC2 

Order 

RD-CC3 

Order 

32 x 32 

7.04E-03 

- 

6.93E-03 

- 

6.16E-04 

- 

8.58E-05 

- 

64 x 64 

2.32E-03 

1.60 

2.20E-03 

1.65 

1.46E-04 

2.08 

6.61E-06 

3.70 

128 x 128 

7.67E-04 

1.60 

6.91E-04 

1.67 

3.82E-05 

1.93 

6.50E-07 

3.35 

256 x 256 

2.54E-04 

1.59 

2.14E-04 

1.69 

1.03E-05 

1.89 

7.62E-08 

3.09 
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Table 8. Iterative convergence of the various hyperbolic RD schemes on irregular grids ( A = 2, C = —1, a = 2, 
b = 1, v = 0.01). Solutions are considered converged when the minimum of ten orders of magnitude reduction 
is achieved for all the nodal-residuals. The Gauss— Seidel relaxation is performed for each Newton iteration 
until at least five orders of magnitude reduction is obtained in the linear residuals. 


Grids 

RD 


RD-JCP2010 

RD- 

-CC2 

RD- 

-CC3 

Newton (Nt) 

GS/Nt 

Newton 

GS/Nt 

Newton 

GS/Nt 

Newton 

GS/Nt 

32 x 32 

3 

68 

3 

48 

3 

72 

10 

76 

64 x 64 

3 

103 

3 

69 

3 

76 

10 

83 

128 x 128 

3 

165 

3 

111 

3 

118 

10 

170 

256 x 256 

3 

273 

3 

190 

3 

249 

10 

361 


rapidly in 10 Newton iterations. Note also that asymptotically the number of linear relaxations increases 
linearly (not quadratically as typically in conventional methods) with the grid sizes. This is a feature of the 
hyperbolic method, arising from the elimination of the second derivatives in the target equation. 5 

To demonstrate the high quality of the predicted solution gradients within the domain by the improved 
RD-CC2 and RD-CC3 schemes, we show results for the 32 x 32 irregular grid shown in Fig. 5b. Accurate 
prediction of the solution gradients within the domain is of great interest in many applications (e.g. turbulent 
flows, reacting flows, etc.). Figure 7 shows the solution gradients predicted by the baseline RD and RD- 
JCP2010 schemes, and compares them with the exact solutions. As can be seen clearly, the predicted results 



(a) u x \ RD (baseline) 


(b) u x ; RD-JCP2010 


(c) u x (exact) 



Figure 7. Comparisons between solution gradients (u x and u y ) predicted by the baseline and RD— JCP2010 
schemes on irregular 32 X 32 grid (see Fig. 5b) with the exact values (A = 2, C = —0.009, a = 2, b = 1, v = 0.01). 


are very oscillatory and inaccurate. On the other hand, as shown in Fig. 8, the gradients predicted by 
the proposed RD-CC2 and RD-CC3 schemes are very accurate and smooth. Note that, our numerical 
experiments suggest that in some cases the noise produced by the baseline RD and RD-JCP2010 schemes 
will not completely disappear in the entire domain even on a highly refined grid. 
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Figure 8. Comparisons between solution gradients (u x and u y ) predicted by the RD- CC2 and RD- CC3 schemes 
on irregular 32 x 32 grid (see Fig. 5b) with the exact values (A = 2, C = —0.009, a = 2, 6=1, v = 0.01). 
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Shown in Figs. 9 and 10 are the reconstructed gradients computed by the quadratic LSQ fits from the 
solution of the RD-CC3 scheme, compared with the gradients predicted directly with the RD-CC3 scheme. 
Clearly the reconstructed gradients are less accurate and deteriorate even more within the domain. This 
comparison emphasizes the difficultly in obtaining accurate solution gradients by reconstruction techniques, 
a common practice with conventional schemes. 



Figure 9. Comparison of u x reconstructed with the quadratic LSQ fit from the solution of the RD— CC3 scheme 
and the one predicted directly with the RD— CC3 scheme on irregular 32 x 32 grid in Fig. 5b. (A = 2, C = —0.009, 
a = 2, 6 = 1, v = 0.01). 


2. Uniform accuracy over Re = ajv 

In this section we demonstrate uniform accuracy predicted by the RD-CC2 and RD-CC3 schemes for both 
the solution and the solution gradients over ranges of Re with different irregular grids. These results are 
shown in Fig. 11, which verifies that the RD-CC2 and RD-CC3 schemes produce, respectively, at a minimum, 
second- and third-order solution for all the variables. Note that in the advection-limit ( v — > 0) the hyperbolic 
advection-diffusion system reduces to simple scalar advection equation, for which the RD-CC2 scheme is 
third-order accurate for all the variables. In the diffusion-limit (y — * oo), the hyperbolic advection-diffusion 
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Figure 10. Comparison of u y reconstructed with the quadratic LSQ fit from the solution of the RD- CC3 
scheme and the one predicted directly with the RD— CC3 scheme on irregular 32 x 32 grid in Fig. 5b. ( A = 2, 
C = -0.009, a = 2, b = 1, v = 0.01). 


24 of 40 


American Institute of Aeronautics and Astronautics 




system becomes hyperbolic diffusion system, for which the RD-CC3 schemes is fourth-order for u but 
remains third-order for the solution gradients. 


RD-CC2, Irregular Grid 



(a) RD-CC2; u 


RD-CC2, Irregular Grid 



RD-CC2, Irregular Grid 



RD-CC3, Irregular Grid 



RD-CC3, Irregular Grid 



RD-CC3, Irregular Grid 



(f) RD-CC3; u y 


Figure 11. Uniform accuracy over ranges of Re with the proposed RD-CC2 and RD- CC3 schemes on irregular 
grids (A = 2, C = —0.009, a = 2, b s= 1, Re = 


Figure 12 shows Newton iteration versus Re for both RD-CC2 and RD-CC3 schemes. It is clear that 
the number of Newton iterations is independent of the Re and grid sizes. We remark that some instabilities 
were observed with the linear relaxation for high-Re cases ( Re > 10 3 ), but found that both RD-CC2 and 
RD-CC3 schemes become stable if the curvature correction term applied to the advection fluxes (i.e. , 6f ) , 
is evaluated with quadratic LSQ reconstructions; that is, u x and u y variables present in the Sf term (see 
Eq. 71) are reconstructed from the solution variable u, instead of being substituted by the gradient variables 
p and q. The jump in the number of Newton iteration shown in Fig. 12 for RD-CC2 is due to this switch, 
which is implemented to avoid instability for high-Re cases. We also remark that if the switch is activated 
for all the Re, the jump will disappear and a similar result shown for RD-CC3 scheme will be recovered. 

A linear convergence, O(N), was also obtained for the ranges of Re. This is shown in Fig. 13. For 
high -Re cases, a better than linear convergence is obtained because of high mesh-Re (i.e., ah/v) values 
used with coarser grids. For example, for the v = 10 -3 case (Re = 1000\/5), linear convergence is obtained 
from the second coarsest grid (mesh-Re < 35). Note that the RD-CC2 and RD-CC3 schemes are efficient 
enough in converging high-Re cases on grids with high mesh-Re, and therefore we did not attempt to fix a 
mesh-Re on these cases. 

3. High aspect-ratio grids 

In practical applications, we are often interested in the gradients at the boundary and with highly anisotropic 
grids, e.g., viscous drag or heat flux predictions in high-Reynolds-number viscous flows. To demonstrate 
the benefit of the hyperbolic schemes for such applications, we consider two highly stretched anisotropic 
irregular grids (see Fig. 14). The y-derivative of u at y = 0 is plotted and compared with the exact solution 
in Fig. 14. The predicted normal gradient is very accurate and smooth on such relatively coarse irregular 
grids. Similar results have been obtained by other schemes, and therefore not shown. 

This exceptional feature of the hyperbolic RD schemes is further illustrated by comparing the normal 
gradients predicted by the baseline RD scheme and a conventional FV scheme on a high aspect-ratio (AR) 

25 of 40 


American Institute of Aeronautics and Astronautics 









RD-CC2, Irregular Grid 


RD-CC3, Irregular Grid 



(a) RD-CC2 


(b) RD-CC3 


Figure 12. Newton iterations used by the proposed RD— CC2 and RD— CC3 schemes on irregular grids over 
ranges of Re and grids (A = 2, C —0.009, a = 2, 6=1, Re = y/b/v). 


RD-CC2, Irregular Grid 



RD-CC3, Irregular Grid 



Figure 13. Average Gauss-Seidel relaxations per Newton iteration used by the proposed RD— CC2 and RD— 
CC3 schemes on irregular grids over ranges of Re and grids (A = 2, C = —0.009, a = 2, 6=1, Re = y/b/u). 
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Derivatives of u at y=0 



(c) 16 x 16 (d) 32 x 32 

Figure 14. Anisotropic, stretched, and irregular grids, and the corresponding predictions of the normal 
gradients along the y = 0 boundary with the baseline hyperbolic RD scheme (A = 2, C = —0.009, a = 2, 6=1, 
v = 0.01). 
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grid. The grid and the results are shown in Fig. 15 where the normal gradients at the y = 0 and the x = 0 
boundaries predicted with the hyperbolic RD scheme are compared with those obtained by the linear LSQ 
reconstruction from the solution of a conventional FV scheme taken from Ref. 8, where the scalar advection- 
diffusion equation is solved by the second-order edge-based advection scheme and the second-order Galerkin 
discretization for diffusion (referred to as Galerkin in the figure). The predicted normal derivative using the 
baseline RD scheme is remarkably smooth and accurate, while the gradients obtained by the conventional 
scheme shows oscillatory behavior and inaccurate predictions. For this problem also, similar results have 
been obtained by other schemes, and therefore not shown. 

These results indicate that even the baseline RD and RD-JCP2010 schemes can produce accurate gra- 
dients on boundaries although they yield oscillations in the interior as shown in the previous section. Also, 
the results demonstrate that the RD-CC2 and RD-CC3 schemes perforin successfully even on irregular 
high -aspect-ratio grids, which are representatives of adapted viscous grids. 



X 


(a) High AR grid (b) u y (c) u x 

Figure 15. Comparisons between the gradients (u x and u v ) predicted with the baseline hyperbolic RD scheme 
and the published result of the conventional Galerkin method 8 on a high aspect-ratio grid with AR = 100; 
the inset figure shows a portion of the grid with 1:1 scaling on both axes. (A = 2, C = 1, a = 1.23, b = 0.12, 
v = 0.1235839795). 


B. Domain with curved geometrical boundaries 


In this section, we verify the accuracy of the RD-CC2 and RD-CC3 schemes for a problem with a curved 
boundary. The problem is a potential flow over a unit circle, taken from Ref. 29 The governing equation is 
the Laplace equation for the stream function, ip: 

dxxi 1 + dyytp = 0 . ( 120 ) 


We solve the above equation by solving the following hyperbolic system of equations: 


d T tp = d x p + d y q, d T p 


1 

Tr 


(d x ip - p ) , 


d T q 


1 

Tr 


(■ dytp - q ) , 


( 121 ) 


where the x— and y-velocity components are related to the gradient of the main variable tp m , i.e. , tp x = — 1 > 
and ip y = u. 

We discretize the hyperbolic potential flow system of equations by the RD-CC2 and RD-CC3 schemes, 
described in Secs. V and VI, and verify their predicted orders of accuracy over a set of eight irregular 
anisotropic triangular grids of 441, 1681, 3721, 6561, 10201, 14641, 25921, and 40401 nodes. The grids are 
composed of linear elements only, and no curved elements are used. The domain and a sample grid (1681 
nodes) is shown in Fig. 16a. We impose the exact solution on the outer domain, and specify ip = 0 at 
the solid boundary nodes to focus on the order of accuracy in the gradients. Note that the gradients at 
boundary nodes are predicted by the numerical schemes. Figure 16 shows the contours of ip x obtained with 
the RD-CC3 scheme on the 1681-node grid. The excellent quality of the predicted derivatives is observed 
in comparison with the exact solution. The predicted streamline is also shown in the same figure. 

The order of accuracy of theses schemes on predicting solution gradients on curved boundaries is also 
verified. The error convergence results for the boundary nodes and the entire domain are shown in Fig. 17. 
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(a) Irregular grid (1681 nodes) 



(c) Predicted streamline on the grid of 1681 nodes. 



(b) Exact (—'i/’x) on the grid of 1681 nodes 



(d) Predicted (—'if’x) on the grid of 1681 nodes. 


Figure 16. Potential flow (ip X x + ^ yy = 0) over a unit circle on irregular and stretched grids with RD- CC3. 
The contours show the x— gradient of the solution variable t/j. 
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Potential flow, Irregular & Stretched Grid Potential flow, Irregular & Stretched Grid 




Figure 17. Potential flow ('ipxx + ^ yy = 0) over a unit circle on several irregular and stretched grids with the 
RD- CC2 and RD- CC3 schemes. 


As it is shown, the design order of accuracy is clearly observed through the boundary nodes. It is also 
evident that the RD-CC3 scheme is more accurate than the RD-CC2 scheme, even on the coarsest grid 
level, indicating that the third-order scheme, which is constructed for linear elements, is, in fact, solving the 
correct problem on the curved boundary that is defined only with linear elements. 


C. Nonlinear advection diffusion problem 

Consider the following nonlinear advection-diffusion equation 

d t u + d x f + d y g = d x {vd x u) + d y {vd y u) + s(ai, y), (122) 

where / = g = u 2 / 2, v = u. The source term s(x, y) is defined by 

s{x, y) = u e (u e x + u e y ) - ( u % ) 2 - ( u e y ) 2 - u e {u e xx + u e yy ), (123) 

^+0), (124) 

where C = —0.009, A = 2.0, = 0.01, and Cg = 1, so that u e is the exact solution to the nonlinear 

advection-diffusion equation (122). Note that Cg must be greater than C in order for the diffusion coefficient 
to be positive. 

We solve the nonlinear advection-diffusion equation in a square domain of [0, 1] x [0, 1], and verify the 
order of accuracy of the solution and solution gradients predicted by the baseline RD, the RD-CC2, and 
the RD-CC3 schemes. The order of accuracy comparison is shown in Fig. 18 along with a more detailed 
information in Tables 9-11, illustrating that the RD-CC2 and RD-CC3 schemes achieve an improved order 
of accuracy for both the solution and the solution gradients on irregular grids. On the other hand, the 
baseline RD scheme suffers from the reduced order of accuracy in the gradients. 


u e = C cos( Airy) exp 


1 - s/1 + 4A 2 tt 2 


2l4l 


IX. Conclusions 

We have presented detailed formulation and implementation procedures for hyperbolic RD schemes for 
general advection-diffusion problems on arbitrary triangular grids. Improved upon the previous work, 5,6 we 
developed new second- and third-order hyperbolic RD schemes that preserve, respectively, quadratic and 
cubic solutions on arbitrary triangular grids. We also developed a higlr-order source term discretization, 
which was found essential in constructing the third-order scheme with linear elements. We showed that the 
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Figure 18. Order of accuracy comparisons between the baseline and the proposed hyperbolic RD schemes on 
irregular grids for the nonlinear advection— diffusion problem (/ = g = u 2 /2, v = u). 


Table 9. L\ error convergence of u for the hyperbolic RD schemes with the SUPG distribution on irregular 
grids, for the nonlinear advection— diffusion problem. 


Grids 

RD 

Order 

RD-CC2 

Order 

RD-CC3 

Order 

32 x 32 

7.93E-05 

- 

8.79E-05 

- 

3.98E-06 

- 

64 x 64 

1.93E-05 

2.04 

2.15E-05 

2.03 

2.43E-07 

4.04 

128 x 128 

4.90E-06 

1.98 

5.46E-06 

1.98 

1.59E-08 

3.93 

256 x 256 

1.22E-06 

2.01 

1.36E-06 

2.01 

1.12E-09 

3.83 


Table 10. L\ error convergence of u x for the hyperbolic RD schemes with the SUPG distribution on irregular 
grids, for the nonlinear advection— diffusion problem. 


Grids 

RD 

Order 

RD-CC2 

Order 

RD-CC3 

Order 

32 x 32 

1.00E-03 

- 

5.66E-04 

- 

5.84E-05 

- 

64 x 64 

3.63E-04 

1.46 

1.59E-04 

1.83 

5.94E-06 

3.30 

128 x 128 

1.44E-04 

1.33 

4.84E-05 

1.72 

7.75E-07 

2.94 

256 x 256 

5.31E-05 

1.44 

1.47E-05 

1.72 

1.12E-07 

2.79 


Table 11. L\ error convergence of u y for the hyperbolic RD schemes with the SUPG distribution on irregular 
grids, for the nonlinear advection— diffusion problem. 


Grids 

RD 

Order 

RD-CC2 

Order 

RD-CC3 

Order 

32 x 32 

1.46E-03 

- 

6.59E-04 

- 

6.17E-05 

- 

64 x 64 

5.62E-04 

1.38 

1.71E-04 

1.95 

6.28E-06 

3.30 

128 x 128 

2.23E-04 

1.33 

4.91E-05 

1.80 

7.94E-07 

2.98 

256 x 256 

8.23E-05 

1.44 

1.44E-05 

1.77 

1.14E-07 

2.80 
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second- and third-order schemes can be constructed by a separate treatment of the advective and diffusive 
terms, enabling development of hyperbolic residual-distribution schemes for the compressible Navier-Stokes 
equations. The superiority of the proposed schemes was further demonstrated with accurate and smooth 
predictions of the gradients within domains, on both flat and curved boundaries. The new second-order 
scheme, which is still defined within a compact stencil, allows an efficient construction of Newton’s method 
with the exact residual Jacobian computed by the automatic differentiation. We demonstrated that the 
resulting implicit solver is capable of reducing the residuals of all equations by at least ten orders of magnitude 
with typically less than 10 Newton iterations for the second-order scheme. For the third-order scheme, we 
showed that the linear solver converges also very rapidly with 10-15 Newton iterations with the second-order 
Jacobian. Because of the hyperbolic reformulation of the advection-diffusion equation, the number of linear 
relaxations per Newton iteration was also shown to increase linearly, not quadratically as for typical diffusion 
problems, with grid refinement. 
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A. Fourier analysis 


In this appendix, we present the Fourier analysis to investigate the stability of relaxation schemes applied 
to the linear system (14). The analysis is performed for the baseline RD scheme only, but the results are 
suggestive and in fact effective to stabilize other schemes. 

Consider a two-dimensional stencil shown in Fig. 19, where we seek the solution at node j, which is 



Figure 19. Schematic of the stencil used in Fourier analysis. 


surrounded by six elements and and six nodes. The grid aspect radio (AR) is defined as h x /h y . We now 
define the vectors U and Q to all the nodes; that is: 


u,= 


Uj 
Vi 

Qj J 


Uz, = 


Uk 

Pk 

Qk 


for k = 1..6, 


(125) 


Q j 


o 

Pj /Tr 
Qj /Tr 


Qfc = 


0 

Pk/T r 

qk/Tr 


for k = 1..6. 


The cell residual 3? T then becomes: 


# T = -^K ?; U. i + Q T 5 T . 


i£T 


For example, for the element containing nodes ( j , 1, 2), we have: 

$T G( yi,2) = _ KjU? _ Kl Ur - K 2 U 2 + ^(Qj + Qr + Q 2 )^yA 


(126) 


(127) 


(128) 
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We remark that K is an element dependent matrix and therefore special care must be taken in evaluating 
the elemental residual $ T . 

Using Eq. (11), we evaluate the residual of the node j, noting that the median dual volume Sj for our 
stencil is simply h x h y . We now start the Fourier analysis by writing the Res_j as 


where C j and C k are: 


6 


; = Cj-u j + ^ 

(129) 

k = 1 


<9Resj 


p j 

- au, ■ 

(130) 

9Res, 


p J 

k dV k ' 

(131) 


We then transform the residual ReSj into the corresponding Fourier modes by replacing the vector Ufc with 
a Fourier mode of the phase angle /? = (/3 X , (3 y ) with j3 € [0, 7r] : 


U r — U i e 


:(/ 3 * 


L +A 


Vk-Vj 

hy 


(132) 


where i = \/— 1, and ( Xj,yj ) and ( Xk,yk ) are the nodal coordinates. Equation 129 can now be written as: 
Res^ = Uj (Cj + e“ i/3 -Ci + e -<(A.+A,)c 2 + e“ i/3 «C 3 + e i/3 *C 4 + e i( ~^ + ^C 5 + e i/3 "C 6 ) (133) 

The solution of Uj using the Jacobi relaxation is therefore become: 

U” +1 = M jb \J] (134) 

where the mass matrix M is defined as 

Mj B = Cj 1 + e -dfe+^)c 2 + e~ iPv C 3 + e^C 4 + e^+^Cs + e i/3 *C 6 ) (135) 


Note that this is equivalent to the Jabobi relaxation applied to the linear system (14) without the right hand 
side, which is irrelevant to the stability analysis. 

The three complex eigenvalues of the matrix, Mjg, are computed for several RD schemes with grid 
AR of 1 and 10 6 . Here we used 50 equally distributed points in /3 X and (3 y directions, for a total of 2500 
points (i.e. , 50 x 50). We first present the results of the first-order Lax-Friedrichs (or Rusanov) scheme 30 
(Fig. 20) for the hyperbolic diffusion problem. The x- and y - axis are the real and imaginary part of the 
matrix eigenvalues. The colors in the plots correspond to the three eigenvalues. This Rusanov scheme is 
stable with GS for any regular grid. Note that the imaginary component of the eigenvalues are very small. 

Figure 21 shows the results for the hyperbolic diffusion system using LDA scheme. Clearly, the results 
exhibit modes with magnitudes larger than one, specially on the real axis. These modes indicate that the 
errors could amplify and may not be damped. Therefore, the LDA scheme is unstable with the Jacobi 
relaxation for our system equation. Similar results are obtained with the Gauss-Seidel relaxation. The 
results are consistent with our developed two-dimensional RD code. We remark that while the figures are 
shown for AR =10° are identical to results obtained on uniform grids with arbitrary AR. Similar results 
are also obtained with a pure hyperbolic problem, the wave equation, u tt = c 2 {u xx +u yy ) (see Fig. 22). The 
results are nearly identical to the hyperbolic diffusion problem shown in Fig. 21 particularly with AR = 1. 
The results indicate that the instability of the LDA scheme is not caused by the first-order hyperbolic system 
method. We have also examined the stability of the LDA scheme for a scalar linear advection equation. Note 
that there is only one eigenvalue for the pure advection problem. These results are presented in Fig. 23, 
which indicate non-damping modes specially for grids with AR = 1. The results confirm that the presented 
instability is not due to the formulation of the diffusion problems using hyperbolic first-order system. 

The SUPG, LW, and a stabilized LDA schemes also exhibit similar instability with the GS relaxation for 
both diffusion and scalar linear advection problems (see Figs. 24-28). Note that the stabilized-LDA scheme 
consists of the LDA scheme and a stabilization term defined as K. 
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(a) Rusanov: AR = 1 (Diffusion) 

Figure 20. Fourier analysis of the Rusanov scheme 



-0.03 - 
-0.04- 

(b) Rusanov: AR = 10 6 (Diffusion) 

for the hyperbolic diffusion problem: v = 1, h x = 0.1. 
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(a) LDA: AR = 1 (Diffusion) 


(b) LDA: AR = 10 6 (Diffusion) 


Figure 21. Fourier analysis of the LDA scheme for the hyperbolic diffusion system: v = 1, h x = 0.1. 
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(a) LDA: Af?, = 1 (Wave) 


(b) LDA: AR = 10 6 (Wave) 


Figure 22. Fourier analysis of the LDA scheme for the hyperbolic wave problem: utt = c 2 {u xx + u yy ), c = 1, 
h x = 0.1. 



0.5 


-*1 


-0.5 


-0.5 


'*** 


0.5 


* * * 5 


(a) LDA: AR = 1 (Advection) 


(b) LDA: AR = 10® (Advection) 


Figure 23. Fourier analysis of the LDA scheme for linear advection equation: a = 2, b = 1, h x = 0.1. 
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(a) SUPG: AR = 1 (Diffusion) 


(b) SUPG: AR = 10 6 (Diffusion) 


Figure 24. Fourier analysis of the SUPG scheme for the hyperbolic diffusion problem: v = 1, h x = 0.1. 



(a) SUPG: Ai? = 1 (Advection) 


(b) SUPG: AR = 10 6 (Advection) 


Figure 25. Fourier analysis of the SUPG scheme for linear advection equation: a = 2, b = 1, h x = 0.1. 
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(a) LW: AR = 1 (Diffusion) 


(b) LW: AR = 10 6 (Diffusion) 


Figure 26. Fourier analysis of the Lax— Wendrof (LW) scheme for the hyperbolic diffusion problem: v = 1 
h x = 0.1. 
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(a) LW: AR = 1 (Advection) 


(b) LW: AR = 10 6 (Advection) 



Figure 27. Fourier analysis of the Lax— Wendrof (LW) scheme for linear advection equation: a = 2, b — 1 
h x = 0.1. 
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(a) Stabilized— LDA: AR = 1 (Diffusion) (b) Stabilized-LDA: AR = 10® (Diffusion) 

Figure 28. Fourier analysis of the stabilized— LDA scheme for the hyperbolic advection— diffusion problem: 
v = 1, h x = 0.1. 


We further examined the instability of these schemes with weighted Jacobi and under-relaxed GS relax- 
ations. The matrix for the weighted Jacobi relaxation, is: 

TsAwjb = (1 — w)I + wMjb, (136) 

where I is an identity matrix and M,;e is defined in Eq. 135. We found that generally these schemes are stable 
with weighted Jacobi or under-relaxed GS relaxations with w = 0.7— 0.8. For instance, Fig. 29 presents the 
stability plots for the linear advection problem using SUPG formulation. The results indicate that the linear 
solver is stable on arbitrary regular grids with cu = 0.7 as the under relaxation value. The same technique 
applied to the hyperbolic diffusion problem as shown in Fig. 30. Similar results were obtained with other 
presented schemes and therefore not repeated here. Note that in the Fourier analysis, effects of the boundary 
points and grid irregularity could not be included. However, we found that the under-relaxation parameter 
of u> = 0.6 — 0.8 is typically sufficient to make the linear relaxation stable on arbitrary anisotropic irregular 
grids for general advection-diffusion problems. 
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